import math
import os
import sys
import dxfgrabber
import simplekml

class Transverse_Mercator(object):   
    def __init__(self,a,f,lat0,lon0,k,FN,FE):
        self.a = a
        self.f = f
        self.lon0 = lon0
        self.lat0 = lat0
        self.k = k
        self.FN = FN
        self.FE = FE

    def proj(self,lat,lon): #横轴墨卡托投影,https://www.iogp.org/bookstore/product/coordinate-conversions-and-transformation-including-formulas/
        e = (2*self.f-self.f**2)**0.5
        n = self.f/(2-self.f)
        B = (self.a/(1+n))*(1+n**2/4+n**4/64)

        h1 = n/2-2/3*n**2+5/16*n**3+41/180*n**4
        h2 = 13/48*n**2-3/5*n**3+557/1440*n**4
        h3 = 61/240*n**3-103/140*n**4
        h4 = 49561/161280*n**4   

        if(self.lat0 == 0):
            Mo = 0
        elif(self.lat0 == math.pi):
            Mo = B*math.pi/2
        elif(self.lat0 == -math.pi):
            Mo = -B*math.pi/2
        else:
            Qo = math.asinh(math.tan(self.lat0))-(e*math.atanh(e*math.sin(self.lat0)))
            βo = math.atan(math.sinh(Qo))
            ξo0 = βo
            ξo1 = h1*math.sin(2*ξo0)
            ξo2 = h2*math.sin(4*ξo0)
            ξo3 = h3*math.sin(6*ξo0)
            ξo4 = h4*math.sin(8*ξo0)
            ξo = ξo0+ξo1+ξo2+ξo3+ξo4
            Mo = B*ξo
        
        Q = math.asinh(math.tan(lat))-(e*math.atanh(e*math.sin(lat)))
        β = math.atan(math.sinh(Q))
            
        η0 = math.atanh(math.cos(β)*math.sin(lon-self.lon0))
        ξ0 = math.asin(math.sin(β)*math.cosh(η0))

        ξ1 = h1*math.sin(2*ξ0)*math.cosh(2*η0)
        ξ2 = h2*math.sin(4*ξ0)*math.cosh(4*η0)
        ξ3 = h3*math.sin(6*ξ0)*math.cosh(6*η0)
        ξ4 = h4*math.sin(8*ξ0)*math.cosh(8*η0)
        ξ = ξ0+ξ1+ξ2+ξ3+ξ4

        η1 = h1*math.cos(2*ξ0)*math.sinh(2*η0)
        η2 = h2*math.cos(4*ξ0)*math.sinh(4*η0)
        η3 = h3*math.cos(6*ξ0)*math.sinh(6*η0)
        η4 = h4*math.cos(8*ξ0)*math.sinh(8*η0)
        η = η0+η1+η2+η3+η4

        self.E = self.FE+self.k*B*η
        self.N = self.FN+self.k*(B*ξ-Mo)

    def r_proj(self,E,N):   #横轴墨卡托投影逆运算,https://www.iogp.org/bookstore/product/coordinate-conversions-and-transformation-including-formulas/
        e = (2*self.f-self.f**2)**0.5
        n = self.f/(2-self.f)
        B = (self.a/(1+n))*(1+n**2/4+n**4/64)

        h1 = n/2-2/3*n**2+5/16*n**3+41/180*n**4
        h2 = 13/48*n**2-3/5*n**3+557/1440*n**4
        h3 = 61/240*n**3-103/140*n**4
        h4 = 49561/161280*n**4   

        if(self.lat0 == 0):
            Mo = 0
        elif(self.lat0 == math.pi):
            Mo = B*math.pi/2
        elif(self.lat0 == -math.pi):
            Mo = -B*math.pi/2
        else:
            Qo = math.asinh(math.tan(self.lat0))-(e*math.atanh(e*math.sin(self.lat0)))
            βo = math.atan(math.sinh(Qo))
            ξo0 = βo
            ξo1 = h1*math.sin(2*ξo0)
            ξo2 = h2*math.sin(4*ξo0)
            ξo3 = h3*math.sin(6*ξo0)
            ξo4 = h4*math.sin(8*ξo0)
            ξo = ξo0+ξo1+ξo2+ξo3+ξo4
            Mo = B*ξo
        
        h1_=n/2-2/3*n**2+37/96*n**3-1/360*n**4
        h2_=1/48*n**2+1/15*n**3-437/1440*n**4
        h3_=17/480*n**3-37/840*n**4
        h4_=4397/161280*n**4

        η_=(E-self.FE)/(B*self.k)
        ξ_=((N-self.FN)+self.k*Mo)/(B*self.k)

        ξ1_=h1_*math.sin(2*ξ_)*math.cosh(2*η_)
        ξ2_=h2_*math.sin(4*ξ_)*math.cosh(4*η_)
        ξ3_=h3_*math.sin(6*ξ_)*math.cosh(6*η_)
        ξ4_=h4_*math.sin(8*ξ_)*math.cosh(8*η_)
        ξ0_=ξ_-(ξ1_+ξ2_+ξ3_+ξ4_)

        η1_=h1_*math.cos(2*ξ_)*math.sinh(2*η_)
        η2_=h2_*math.cos(4*ξ_)*math.sinh(4*η_)
        η3_=h3_*math.cos(6*ξ_)*math.sinh(6*η_)
        η4_=h4_*math.cos(8*ξ_)*math.sinh(8*η_)
        η0_=η_-(η1_+η2_+η3_+η4_)

        β_=math.asin(math.sin(ξ0_)/math.cosh(η0_))
        Q_=math.asinh(math.tan(β_))

        temp_Q=Q_
        while(1):
            old_Q=temp_Q
            temp_Q=Q_+e*math.atanh(e*math.tanh(temp_Q))
            if(math.fabs(old_Q-temp_Q)<10**(-10)):
                break
        self.lat=math.atan(math.sinh(temp_Q))
        self.lon=self.lon0+math.asin(math.tanh(η0_)/math.cos(β_))
  
    def geo2xyz(self,lat,lon,h):
        ee = 2*self.f-self.f**2
        v = self.a/(1-ee*math.sin(lat)**2)**0.5

        self.x = (v+h)*math.cos(lat)*math.cos(lon)
        self.y = (v+h)*math.cos(lat)*math.sin(lon)
        self.z = ((1-ee)*v+h)*math.sin(lat)
  
    def xyz2geo(self,x,y,z):
        ee = 2*self.f-self.f**2
        
        self.lat = 0
        temp_l = (x**2+y**2)**0.5
        while(1):
            v = self.a/(1-ee*math.sin(self.lat)**2)**0.5
            temp_lat = math.atan2(z+ee*v*math.sin(self.lat),temp_l)
            if(math.fabs(temp_lat-self.lat)<10**(-10)):
                break
            self.lat = temp_lat
        self.lon = math.atan2(y,x)
        self.H   = x/math.cos(self.lat)/math.cos(self.lon)-v

    def test(self):
        print("经度=",math.degrees(self.lon))
        print("纬度=",math.degrees(self.lat))
        print("高度=",self.H)
        #print("x=",self.x)
        #print("y=",self.y)
        #print("z=",self.z)

class wgs84(Transverse_Mercator):
    def __init__(self):
        self.a = 6378137
        self.f = 1/298.2572236

class sld99(Transverse_Mercator):
    def __init__(self):
        self.a = 6377276.345
        self.f = 1/300.801699999587
        self.lon0 = 1.40973233578035
        self.lat0 = 0.12218127735184
        self.k = 0.9999238418
        self.FN = 500000
        self.FE = 500000

        self.tx = -0.2933
        self.ty = 766.9499
        self.tz = 87.7131
        self.rx = 0.000000948799766478606           #0.1957040"
        self.ry = 0.000008217920113668750           #1.6950677"
        self.rz = 0.000016837657199936800           #3.4730161"
        self.dS = -0.000000039338
    
    def to_wgs84(self,E,N,h=0): 
        point_a = wgs84()
        self.r_proj(E,N)
        self.geo2xyz(self.lat,self.lon,h) 

        M = 1+self.dS
        x = M*(+1       *self.x  -self.rz *self.y +self.ry  *self.z)+self.tx
        y = M*(+self.rz *self.x  +1       *self.y -self.rx  *self.z)+self.ty
        z = M*(-self.ry *self.x  +self.rx *self.y +1        *self.z)+self.tz
        
        point_a.xyz2geo(x,y,z)

        return point_a

def file_name(file_dir,file_type):
    L=[]
    for root,dirs,files in os.walk(file_dir):
        for file in files:
            if os.path.splitext(file)[1] == file_type:
                L.append(os.path.join(root, file))
    print(L)
    return L

def dxf2kml(file_name):
    dwg = dxfgrabber.readfile(file_name)
    kml = simplekml.Kml()
    P = sld99()

    style = simplekml.Style()
    kml.AltitudeMode = 'relativeToGround'

    nn = 128 #以直代弧分段数

    for e in dwg.entities:
        
        if e.dxftype == 'LINE':
            LL = []
            LL.append(e.start)
            LL.append(e.end)
            for i in range(len(LL)):
                Temp = P.to_wgs84(LL[i][0],LL[i][1])
                LL[i] = (math.degrees(Temp.lon),math.degrees(Temp.lat))
            lin = kml.newlinestring(name = "LINE",coords = LL)
            lin.style = style
            lin.style.linestyle.color = e.true_color[2:] if e.true_color  else 'FFFFFFFF'
        if e.dxftype == 'MTEXT':
            Temp = P.to_wgs84(e.insert[0],e.insert[1])
            pnt = kml.newpoint()
            pnt.name = e.raw_text
            pnt.coords = [(     math.degrees(Temp.lon),math.degrees(Temp.lat)   )]
        if e.dxftype == 'CIRCLE':
            LL = []
            for i in range(0,nn):
                LL.append((e.center[0]+e.radius*math.cos(math.radians(i*math.pi*2/nn)) , e.center[1]+e.radius*math.sin(math.radians(i*math.pi*2/nn)) ))
            for i in range(len(LL)):
                Temp = P.to_wgs84(LL[i][0],LL[i][1])
                LL[i] = (math.degrees(Temp.lon),math.degrees(Temp.lat))
            lin = kml.newlinestring(name = "LINE",coords = LL)
            lin.style = style
            lin.style.linestyle.color = e.true_color[2:] if e.true_color  else 'FFFFFFFF'            

        if e.dxftype == 'ARC':
            LL = []
            angle = e.end_angle - e.start_angle
            for i in range(0,nn):
                LL.append(( e.center[0]+e.radius*math.cos(math.radians(e.start_angle+i*angle/nn)) , e.center[1]+e.radius*math.sin(math.radians(e.start_angle+i*angle/nn)) ))
            for i in range(len(LL)):
                Temp = P.to_wgs84(LL[i][0],LL[i][1])
                LL[i] = (math.degrees(Temp.lon),math.degrees(Temp.lat))
            lin = kml.newlinestring(name = "LINE",coords = LL)
            lin.style = style
            lin.style.linestyle.color = e.true_color[2:] if e.true_color  else 'FFFFFFFF'

        if e.dxftype == 'LWPOLYLINE':
            
            LL = []
            
            for m in range(0,e.__len__()-1):
                if e.bulge[m] == 0:
                    LL.append(e.points[m])
                else:
                    l = ((e.points[m][0]-e.points[m+1][0])**2+(e.points[m][1]-e.points[m+1][1])**2)**0.5
                    h = e.bulge[m]*l/2
                    R = h/2+l**2/8/h

                    sign = 1 if e.bulge[m] < 0 else -1
                    dsign = 1 if math.fabs(e.bulge[m]) >1 else -1
                    d = math.fabs(R-h)
                    cx = (e.points[m][0]+e.points[m+1][0])/2 + d*math.cos(math.atan2((e.points[m+1][1]-e.points[m][1]),(e.points[m+1][0]-e.points[m][0]))+dsign*sign*math.pi/2)
                    cy = (e.points[m][1]+e.points[m+1][1])/2 + d*math.sin(math.atan2((e.points[m+1][1]-e.points[m][1]),(e.points[m+1][0]-e.points[m][0]))+dsign*sign*math.pi/2)

                    starta = math.atan2(e.points[m][1]-cy,e.points[m][0]-cx)
                    enda = math.atan2(e.points[m+1][1]-cy,e.points[m+1][0]-cx)

                    angle = enda - starta
                    
                    if angle < 0:
                        angle = angle+math.pi*2
                    
                    if dsign == 1:
                        if angle < math.pi:
                            angle = math.pi*2-angle
                    else:
                        if angle > math.pi:
                            angle = math.pi*2-angle   

                    for i in range(0,nn):
                        LL.append((cx-sign*R*math.cos(starta-sign*i*angle/nn),cy-sign*R*math.sin(starta-sign*i*angle/nn)))

            if e.is_closed:
                LL.append(e.points[0])
        
            for i in range(len(LL)):
                Temp = P.to_wgs84(LL[i][0],LL[i][1])
                LL[i] = (math.degrees(Temp.lon),math.degrees(Temp.lat))
            lin = kml.newlinestring(name = "LINE",coords = LL)
            lin.style = style
            lin.style.linestyle.color = e.true_color[2:] if e.true_color  else 'FFFFFFFF'

    kml.save(file_name[:-4]+".kml")

def main():
    for fname in file_name(os.path.split(os.path.realpath(__file__))[0],".dxf"):
    #for fname in file_name(os.path.dirname(sys.executable),".dxf"):
        dxf2kml(fname)
        
    
if __name__ == '__main__':
    main()        
    